1 Summary

Examine biodiversity risk by taxonomic group.

2 Data sources

• IUCN species API: IUCN. (2019). The IUCN Red List of Threatened Species. Version 2019-2.
• IUCN species shapefiles: IUCN. (2019). The IUCN Red List of Threatened Species. Version 2019-2. Retrieved August 2019, from http://www.iucnredlist.org
• BirdLife International shapefiles: BirdLife International and Handbook of the Birds of the World. (2018). Bird species distribution maps of the world. Version 7.0. Available at http://datazone.birdlife.org/species/requestdis
• World Database on Protected Areas: IUCN, & UNEP-WCMC. (2018, June). The World Database on Protected Areas (WDPA). Retrieved June 9, 2018, from Cambridge, UK: UNEP-WCMC website: www.protectedplanet.net
• Marine Ecoregions of the World: Spalding, M. D., Fox, H. E., Allen, G. R., Davidson, N., Ferdaña, Z. A., Finlayson, M. A. X., … others. (2007). Marine ecoregions of the world: A bioregionalization of coastal and shelf areas. BioScience, 57(7), 573–583.

3 Methods

spp_maps <- read_csv(file.path(dir_data, 
                               sprintf('spp_marine_maps_%s.csv', api_version)),
                     col_types = 'ddcicccl') %>%
  filter(comp_assessed == TRUE) %>%
  mutate(dbf_file = basename(dbf_file)) %>%
  select(iucn_sid, dbf_file)
spp_info <- read_csv(file.path(dir_o_anx, 
                               sprintf('iucn/spp_info_from_api_%s.csv', api_version))) %>%
  select(-sciname)

taxa_list <- spp_maps %>%
  left_join(spp_info, by = 'iucn_sid') %>%
  select(iucn_sid, kingdom, phylum, class, order, family, dbf_file)
fam_list <- taxa_list %>%
  group_by(phylum, class, order, family, dbf_file) %>%
  summarize(n_spp = n())
order_list <- taxa_list %>%
  group_by(phylum, class, order, dbf_file) %>%
  summarize(n_spp = n())
class_list <- taxa_list %>%
  group_by(phylum, class, dbf_file) %>%
  summarize(n_spp = n())

shp_to_taxa <- read_csv(file.path(dir_setup, 'raw', 'shps_to_taxa.csv'))

taxa_sums <- data.frame(
  unweighted = list.files(file.path(dir_o_anx, 'taxa_sums_mol_2019'),
                          pattern = sprintf('cell_sum_comp_%s.csv', api_version),
                          full.names = TRUE),
  rr_weighted = list.files(file.path(dir_o_anx, 'taxa_sums_mol_2019'),
                           pattern = sprintf('cell_sum_comp_rrweight_%s.csv', 
                                             api_version),
                           full.names = TRUE)) %>%
  mutate(shp_base = basename(unweighted),
         shp_base = str_replace_all(shp_base, '_cell_sum.+|_part.+', '')) %>%
  left_join(shp_to_taxa, by = 'shp_base')

# write_csv(taxa_sums %>% select(shp_base, taxon) %>% distinct(),
#           file.path(dir_setup, 'raw', 'shps_to_taxa_raw.csv'))

spp_maps <- read_csv(file.path(dir_data, sprintf('spp_marine_maps_%s.csv', api_version)),
                     col_types = 'ddcicccl') %>%
  filter(comp_assessed == TRUE) %>%
  mutate(shp_base = tolower(dbf_file) %>%
           basename() %>%
           str_replace_all(., '_part.+|.dbf', '')) %>%
  select(iucn_sid, shp_base) %>%
  distinct() %>%
  group_by(shp_base) %>%
  summarize(n_spp = n()) %>%
  left_join(taxa_sums, by = 'shp_base') %>%
  select(taxon, n_spp) %>%
  distinct() %>%
  group_by(taxon) %>%
  summarize(n_spp = sum(n_spp)) %>%
  arrange(desc(n_spp))

reload <- TRUE ### reload applies to building rasters, not the plot

rast_base <- raster::raster(file.path(dir_spatial, 'cell_id_mol.tif'))

land_poly <- sf::read_sf(file.path(dir_spatial, 'ne_10m_land', 
                                   'ne_10m_land_no_casp.shp')) %>%
  st_transform(crs(rast_base))

taxa_gps <- c("Cnidaria",
              "Tracheophyta",
              "Chordata: Actinopterygii",
              "Chordata: Mammalia",
              "Chordata: Chondrichthyes",
              "Chordata: Aves",
              "Chordata: Reptilia",
              "Echinodermata")

map_list <- vector('list', length = length(taxa_gps)) %>%
  setNames(taxa_gps)

for(i in seq_along(taxa_gps)) { ### i <- 6
    
  taxon_gp <- taxa_gps[i]
  taxon_txt <- taxon_gp %>% 
    tolower() %>% 
    str_replace_all('chordata: |[^a-z]+', '')

  taxa_info <- taxa_sums %>%
    filter(taxon == taxon_gp)
  
  taxa_rast_file <- file.path(dir_o_anx, 'taxa_rasts_mol_2019', 
                              paste0(tolower(taxon_txt), '_rr_comp.tif'))
  
  if(!file.exists(taxa_rast_file) | reload == TRUE) {
  
    cat_msg('Building: ', taxa_rast_file)

    rr_weighted_files <- taxa_info$rr_weighted %>% unique()
    
    taxon_df_rr <- parallel::mclapply(rr_weighted_files, mc.cores = 12,
        FUN = function (x) { ### x <- rr_weighted_files[1]
          read_csv(x, col_types = 'ddd___d___')
        }) %>% 
      bind_rows() %>%
      group_by(cell_id) %>%
      summarize(mean_risk_sum = 1/sum(sr_rr_risk) * sum(mean_risk * sr_rr_risk),
                pct_threat_sum = sum(sr_rr_threatened) / sum(sr_rr_risk)) %>%
      select(cell_id, mean_risk_sum, pct_threat_sum)
    
    mean_rr_rast <- raster::subs(x = rast_base, y = taxon_df_rr,
                                 by = 'cell_id', which = 'mean_risk_sum')
    
    raster::writeRaster(mean_rr_rast, filename = taxa_rast_file, overwrite = TRUE)
    
  } else {
    cat_msg('File exists: ', taxa_rast_file)
    mean_rr_rast <- raster(taxa_rast_file)
  }

  mean_rr_df <- mean_rr_rast %>% 
    aggregate(fact = 2) %>%
      ### should look OK for tiny maps?
    rasterToPoints() %>% 
    as.data.frame() %>%
    setNames(c('long', 'lat', 'value'))

  cat_msg('Building plot for ', taxon_gp)
  ### using geom_tile so I can add a size to pixels, thought it is
  ### significantly slower...
  x <- ggplot(mean_rr_df) +
    geom_tile(aes(long, lat, fill = value, color = value), 
              size = .5,
              show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map(base_size = 9) +
    theme(plot.margin = unit(c(.05, .10, .05, .05), units = 'cm'),
          panel.background = element_rect(fill = 'grey80', color = NA),
          axis.title = element_text(face = 'bold', hjust = 1, 
                                    vjust = 0, size = 7,
                                    margin = margin(t = 0, r = 0, 
                                                    b = 0, l = .1, unit = 'cm')),
          axis.title.x = element_blank()) +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = risk_cols, 
                         values = risk_vals, 
                         limits = c(0, 1),
                         labels = risk_lbls, 
                         breaks = risk_brks) +
    scale_color_gradientn(colors = risk_cols, 
                         values = risk_vals, 
                         limits = c(0, 1),
                         labels = risk_lbls, 
                         breaks = risk_brks) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0)) +
    labs(y = taxon_gp)

  map_list[[i]] <- x
  
}

  ### set up a dataframe of values to craft a color bar using geom_segment
  colorbar_df <- data.frame(x = seq(0, 1, .001), y = -1)

  colorbar <- ggplot(colorbar_df) +
    ggtheme_plot(base_size = 7) +
    geom_segment(aes(x = x, xend = x, color = x), 
                 y = 0, yend = 1, size = 1,
                 show.legend = FALSE) +
    scale_color_gradientn(colors = risk_cols, values = risk_vals, 
                          limits = c(0, 1),
                          labels = risk_lbls, breaks = risk_brks) +
    theme(axis.text.y = element_blank(), 
          axis.title  = element_blank(),
          panel.grid.major.y = element_blank(),
          plot.margin = unit(c(.05, 0, .05, 0), units = 'cm')) +
    scale_x_continuous(labels = risk_lbls, breaks = risk_brks, limits = c(0, 1),
                       expand = expand_scale(.25)) +
    scale_y_continuous(limits = c(0, 1))

cat_msg('done, putting together the cowplot')

y <- cowplot::plot_grid(plotlist = map_list, align = 'hv', ncol = 2)

z <- cowplot::plot_grid(y, colorbar, 
                        rel_heights = c(20, 1),
                        align = 'v', ncol = 1)
plot_w <- 4.75
plot_h <- 6.1

taxa_map_file <- file.path(dir_git, 'ms_figures/fig2_taxa_maps_mol.tiff')
taxa_map_file <- file.path(dir_git, 'ms_figures/fig2_taxa_maps_mol.png')

ggsave(plot = z,
       filename = taxa_map_file,
       width = plot_w, height = plot_h, units = 'in', dpi = 300)

spp_ranges <- read_csv(file.path(dir_data, 
                                 sprintf('iucn_spp_range_area_%s_mol.csv', api_version)),
                       col_types = 'dd__') %>%
  distinct()

taxa_ranges <- read_csv(file.path(dir_data, 
                                  sprintf('spp_marine_maps_%s.csv', api_version)),
                        col_types = 'ddcicccl') %>%
  filter(comp_assessed) %>%
  mutate(shp_base = basename(tolower(dbf_file)),
         shp_base = str_replace_all(shp_base, '.dbf|_part.+', '')) %>%
  left_join(spp_ranges, by = 'iucn_sid') %>%
  left_join(shp_to_taxa, by = 'shp_base') %>%
  filter(comp_assessed &!is.na(taxon)) %>%
  mutate(log_range = log10(range_km2))

x <- ggplot(taxa_ranges, aes(x = log_range, ..scaled..)) +
  ggtheme_plot() +
  theme(panel.grid.major = element_blank(),
        axis.text.y      = element_blank(),
        axis.title.y     = element_blank(),
        strip.text.y = element_text(angle = 0, hjust = 0)) +
  geom_density(color = 'grey20', size = .25, fill = 'grey30') +
  scale_x_continuous(limits = c(0, NA),
                     breaks = c(0:8),
                     labels = 10^(0:8)) +
  labs(x = 'Range, km^2') +
  facet_grid(taxon ~ .)

# ggsave(filename = file.path(dir_git, 'ms_figures', 'fig_SI_taxa_ranges.png'), 
#            width = 4.75, height = 3, units = 'in', dpi = 300)

print(x)